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Abstract. We use a dynamical systems approach to investigate Bianchi type I and 
II universes in quadratic theories of gravity. Due to the complicated nature of the 
equations of motion we focus on the stability of exact solutions and find that there exists 
an isotropic FRW universe acting as a past attractor. This may indicate that there is 
an isotropisation mechanism at early times for these kind of theories. We also discuss 
the Kasner universes, elucidate the associated centre manifold structure, and show that 
there exists a set of non-zero measure which has the Kasner solutions as a past attractor. 
Regarding the late-time behaviour, the stability shows a dependence of the parameters 
of the theory. We give the conditions under which the de Sitter solution is stable and 
also show that for certain values of the parameters there is a possible late-time behaviour 
with phantom-like behaviour. New types of anisotropic inflationary behaviour are found 
which do not have counterparts in general relativity. 



1. Introduction 

In this paper we are going to study the dynamical evolution of a class of anisotropic 
universes in quadratic theories of gravity. These extensions of general relativity (GR) pro- 
vide guidance as to the possible effects of quantum corrections to the Einstein equations. 
They permit us to investigate the possible effects on singularity formation, inflation, and 
the expansion dynamics of the early universe. Past studies of these extensions have focused 
on the isotropic Friedmann metrics, where it is sufficient to consider only the effects of an 
R 2 term in the gravitational lagrangian to the field equations pQ. However, the R 2 con- 
tribution has fairly predictable cosmological consequences because the resulting quadratic 
vacuum theory is conformally equivalent to GR with a scalar field moving in a potential 
with a single asymmetric minimum [2j [3] . This type of solution has been well studied in 
connection with inflation 4J and in the situation of pure power-law lagrangians of the form 
R n , [3 [7l [8] . The addition of the R I1V R ,J ' V Ricci term to the lagrangian in the case of 
an anisotropic universe creates a much richer diversity of cosmological behaviours that are 
harder to summarise in terms of simple modifications of the general-relativistic situation. 
The effective stresses that are contributed to the field equations by the quadratic Ricci terms 
in the lagrangian can mimic a wide range of fluids which violate the strong and weak energy 
conditions. This allows completely different behaviour to occur than is found in GR or its 
quadratic extension with pure R 2 contributions. In particular, we have shown elsewhere 
that the cosmic no-hair theorems no longer hold: vacuum universes with positive cosmo- 
logical constant do not necessarily approach de Sitter, but can inflate anisotropically [9]. 
Moreover, the addition of quadratic Ricci terms to the lagrangian can create cosmological 
models which have no counterpart in the GR limit of the theory. 

In what follows, we are going to widen this investigation by considering the global dy- 
namics of the Bianchi type I and II universes in the presence of quadratic Ricci and scalar 
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curvature contributions to the lagrangian. Whilst these anisotropic universes contain some 
mathematical simplifications, they include spatially homogeneous expanding universes with 
isotropic (type I) and anisotropic (type II) spatial 3-curvatures as well as expansion shear 
anisotropy. To this end, we consider the gravitational action 

(1) S G = ^- I d i x^/\i\(R + aB 2 +l3R^R^-2A). 

Variation of this action leads to the following generalised Einstein equations (see, e.g., [TP]): 

(2) + + Ag^ = kT„ v , 

where T^„ is the energy-momentum tensor of the ordinary matter sources, which in this 
paper we will assume to be zero, for simplicity, 

1 

2 J 



(3) G M „ = R^ v — —Rg^y, 



$ M „ = 2aR ( R„ v - -Rg^ ) + (2a + (3) (g^D - V M V„) R 



(4) +0D Ut MV - ^Rg^j + 2(3 \R^ vp - \g^R*}j R ap , 

with □ = V M V M and g^ v is the metric tensor and A the cosmological constant. The effective 
stress tensor incorporates the deviation from regular Einstein gravity introduced by the 
quadratic terms in the action, and we see that a = (3 = implies <I> M „ = 0, although the 
converse need not be true. 

Similarly, by rescaling (a, (3) i— > n(ot,(3) we can consider the limit k — ► oo. This limit, 
whenever it exists, corresponds to the case where the action only contains quadratic terms, 
hence eq. reduces to (for vacuum, = 0) 

(5) $ M „ - 0. 

It is important to consider the vacuum solutions to the pure quadratic theories; that is, 
solutions to eq.©. First, we note that any Einstein metric with R^ v = Xg^, necessarily 
obeys = 0. In particular, this implies that any GR vacuum solution (G^ v = 0) will 
also be a vacuum solution to the quadratic theory (G^v + = 0). This means that 
there will be chaotic dynamics in the quadratic theory because it is present in a number 
of general-relativistic vacuum Bianchi models (types Vl^jyg, VIII, and IX) near an initial 
curvature singularity; however it might not be stable in solutions to the quadratic theories 
in the way that it appears to be for these solutions in GR. Second, there are some specific 
solutions to <!> M „ = 0, which are not solutions to G^ u = 0, and so have no counterparts 
in the (a,/?) — ► (0,0) general-relativity limit of the quadratic theory. This feature plays 
an important role in the evolution of universes that arise in these theories. For example, 
one such intrinsically quadratic solution is a Friedmann-Robertson- Walker (FRW) universe 
which has expansion dynamics similar to those of a radiation-dominated universe in GR. 
This means that if we take a FRW universe with spatial curvature parameter fc, there is 
a simple solution to — given by the metric that solves the field equations of GR 
(G^ = kT^v) in the presence of black-body radiation: 



(6) ds 2 = -At 2 + (t- kt 2 ) 



1 — fcr 2 



2 (d<^ + shV 4>A6 Z ) 



Note that despite the fact that = this quadratic universe behaves as if it is radiation- 
dominated. This radiation-like universe seems to play a special role in the theory as it will 
be found to act as a past attractor for the anisotropic universes we study. Other interesting 
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conclusions can also be deduced from this simple solution. We see, for example, that unlike 
for the GR case, in the quadratic theory there is a closed (k = 1) FRW vacuum solution 
which recollapses after a finite time (t = 1) and a flat (k = 0) vacuum solution. 

There are some further geometrical observations regarding solutions to these quadratic 
theories that will be useful. In 4D spacetimes we can use the Weyl curvature invariant and 
the Euler density, E, defined by, 

n rt^vpa p Wvpo nn r>liv . \ t>2 

(7) E = R^R^-AR^R^+R 2 , 

to replace the quadratic Ricci invariant in the action, since 

aR 2 + PR^ST = J (3a + (3)R 2 + § (C^ pa C^ pa - E) . 

Since integration over the Euler density is a topological invariant, the variation of E will 
not contribute to the equations of motion. From this we see that there are two special cases 

a. (3a + 0) = : the conformally invariant Bach- Weyl gravity theory [11] . 

b. (3 = : A special case of f(R) gravity. 

Since the equations of motion change their structure in these two special cases, a separate 
analysis is required for each. So, in what follows we will assume that 

3a + /3^0, (3^0. 

On the other hand, our formulation of the equations of motion is well defined in the k — > oo 
limit (as explained above). 

In our analysis, we will adopt the dynamical systems approach by introducing expansion- 
normalised variables. This approach has proven to be extremely successful in the analysis of 
spatially homogeneous Bianchi universes in GR [12L 113] . However, by considering quadratic 
theories of gravity there are some extra complications added to the formalism, as the universe 
can bounce or recollapse at expansion minima and maxima, which leads to infinities in 
the expansion-normalised variables. Here, we will discuss these infinities and explain how 
a sequence of expanding and contracting phases can be considered. On the other hand, 
initial singularities and future asymptotes are all finite in these variables, and the dynamical 
systems approach is ideal for studying such regimes. 

In a recent paper by Leach et al [14] , a dynamical systems approach to the local rotational 
symmetric Bianchi type I universes was considered within a class of f(R) theories. In 
particular, they were interested in the shear dynamics of such universes, and, interestingly, 
they found that there is a possibility for an isotropic singularity. The introduction of a 
RfivR^" term introduces extra shear degrees of freedom and their higher time derivatives, 
which considerately complicate the equations of motion. However, in spite of the fact that 
the case we study is is a quite separate, and we consider more general Bianchi models, we 
also find the possibility for an isotropic singularity, as did Cotsakis et al [15] in their study 
of Bianchi type IX universes in quadratic theories. Indeed, we will argue that this isotropic 
singularity is past stable for all Bianchi models. 

2. Equations of motion 
Our starting point is the generalised vacuum field equations: 

G^ v + <5>^ + K 9pu = 0. 
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We note that the GR-limit can be obtained by letting (a, (3) — > (0,0). Now consider the 
spatially homogeneous Bianchi metrics. We can always write their metric line-elements as 

ds 2 = -dt 2 + 5 ab u> a u> b , 

where u) a is a triad of one-forms obeying 

(d W °) x = -~C bc u; 6 Au; c , 

the C a bc depend only on time and at any moment of constant time are the structure con- 
stants of the Bianchi group-type under consideration, and (— )j_ means projection onto the 
spatially-homogeneous hypersurfaces. The structure constants can be split into a vector 
part a b and a symmetric tensor n ab in the standard way by 

C\ c = e bcd n da - S a b a c + S a c a b . 

Using the Jacobi identity, a b is in the kernel of n ab 

n ab a b = 0. 

The different structures of n ab and a b define the various Bianchi types (see, e.g., [T^]). 

Defining the time-like hypersurface-orthogonal vector u = d/dt, we can define the Hubble 
scalar, if, and the shear, a ab , as follows: 

1 

3' 



H = CT a fc = U( a; 6) _ H5 a b 



We will also restrict attention to cosmological models where the shear is diagonal, so we can 
write 

a ab = diag(-2er + , er + + V3(t_,(t + - V3a-). 

We are interested in the Bianchi type I and type II models for which a b — and we can 
write 

n a b = (type I), n ab = diag(mi, 0, 0) (type II). 

Wc define the dimcnsionlcss expansion-normalised variables by scaling out appropriate 
powers of H 

R 1 P 

B = 77, , m,,, , X 



Q=^ ^2 = 4- Ov = r4>, N= 7111 



(8) S± = ^, E±i = S± 2 = 

Note the presence of time derivatives of the variables Q2 and £+2; this reflects the 4*' l -order 
time derivatives in the field equations of the quadratic theorjo; \ IS a constant. We also 
introduce the dynamical time variable r by 

and we assume that the cosmological constant is positive: tt\ > 0. 



^The variable Q2 can be related to the statefinders g and j (see e.g. 1231 ). by Q2 = j + 3q + 2. 
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The equations of motion and the constraint ('Friedmann-likc') equation are: 



(9) B' = -2QB, 

= -2Qfl A , 

= -(Q + 1 + 4E+)JV, 

= -2Q 2 + Q 2 , 

= -3(Q + 2)Q 2 -^(Q + 2)Q-^£^l + £ 2 -r! A - 

-|(1 + 2 X )£ 4 - 1(8 + *)£? - (4 - X )(S • Si) 
-1(4 - x)(3E 2 + 2E • £ 2 + 2QE 2 ) - (1 + 2Q)iV 2 



(10) n' A 

(11) JV' 

(12) Q' 

Q' 2 



-N- 



(13) 

(14) £' ± 

(15) E' ±1 



+N 2 



(1 + 8 X )^V 2 + 5(13 + 3 X )E 2 + 8(2E+ - E +1 ) + (1 - X )E 2 



— QE ± + E±i, 
-2QE±i + E±2, 

-3(Q + 2)E +2 4 



E+i 

X 



[B - (11 X - 8) + 4Q(1 - X ) + 4E 2 (1 + 2 X )] 



(16) 



+ ^ [3S + (4 - x)(6 + Q 2 + 7Q) + 4(1 + 2 X )(3E 2 + 2E • £i)l 

X 

--N 2 \B + 8 + 4Q - 4(1 + Sx)^ 2 ] 
X 

~-N 2 [(1 + 15x)(E+ + E +1 - 4E 2 ,) + 4(1 - X )E 2 _] 
E'_ 2 = -3(Q + 2)E_ 2 + ^i [B-(llx-8) + 4Q(l-x)+4E 2 (l + 2x) 

+ — [3B + (4 - x)(6 + Q 2 + 7Q) + 4(1 + 2 X )(3E 2 + 2E • Ei)] 
X 

_ 4 (!-x) JV 2/ E _ + s _ i _ 8S _ S+ ). 

X 

These equations are subject to the constraint: 

= B(l-tt A -E 2 -7V 2 ) + 12Q-2Q 2 + 4Q 2 -(4-x)(3 + 2Q)E 2 
-6(1 + 2 X )E 4 - x(E 2 - 2S • E 2 ) + 4(2 + *)(£■ Si) 



(17) 



(18) 



+47V 2 



-(1 + 8 X )N 2 + 1 + (1 + 15x)E 2 h + 8E+ + (1 - x)E 2 



where, we have introduced the short-hand notation E„ = (£ + „,E_„) and (E„ • E m ) = 
£-t-n£+m -(- E_ n S_ m . 

There are two points to note. First, the parameter Q is related to the usual deceleration 
parameter q via 

q = -(l + Q). 

Second, the variable B measures how greatly the quadratic part of the lagrangian dominates 
over the general-relativistic Einstein-Hilbert term R— 2 A. In particular, the larger the value 
of B, the "closer" we are to GR. The B = case corresponds to a purely quadratic lagrangian 
theory whose equations of motion reduce to = 0. 
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The equations of motion define a dynamical flow in a large phase volume, and the be- 
haviour can be analysed qualitatively by standard techniques from the theory of ordinary 
differential equations. Of particular interest are the exact solutions which define the critical 
points of the system and their stability at small and large times. 

3. Solutions and their stability 

3.1. The de Sitter solution: dS. The de Sitter solution is characterised by the critical 
points where 

Q = Q-2 = S± = S±i = S ±2 = N = 0, Ov = 1, 0. 
Its stability is determined by the eigenvalues 



„,- 1 ,-^ 1± ; 1 _-j,- 3[x21 ,-|( 1 ^ 1+ i!^i!) M 

The zero eigenvalue corresponds to the variable B and appears because these solutions are 
a one-parameter family. The last three eigenvalues arise from the shear equations and come 
therefore in pairs. We also see that the limit \ — > is ill-defined for these eigenvalues and 
these arise therefore only for theories with a non-zero 
We note that 

B>0^(3a + P)>0, * + 2(4- X ) <0 ^l + 2A(4q + /?) < q _ 

X P 
In particular, this means that whenever a — 0, [3 > the de Sitter solution has two unstable 
modes. These modes stem from the shear equations and are therefore not present for the 
FRW models studied earlier. On the other hand, if (3 < and (3a + (3) > then the de 
Sitter solution is stable. 

A related set of equilibrium points is 

Q = Q 2 = S± = £±i = S± 2 = N = B = 0, n A = constant ^ 0. 

The eigenvalues for this line bifurcation is the same as for the de Sitter solution restricted 
to B = 0. Consequently, there exist two zero eigenvalues for this set of solutions. However, 
it can be shown that these solutions are always unstable. 

3.2. The Kasner circle: K.. The Kasner circle of equilibrium points is given by 
Q = —3. (S + , £_) = (cos^>,sin0), (S+i, £-i) = — 3(cos </>, sin </>), 
(19) (£ +2 ,£_ 2 ) = 18(cos0,sin^), Q A = B = N = 0. 

The exact Kasner solutions are completed by integrating the B equation so that B = 
B e 6r ; however, here we will mostly be interested in the limit r — > — oo for which case the 
equilibrium points suffice. 

Stability is determined by the eigenvalues: 

0[x2], 2(1 -2 cos(/>), 6[x7]. 



2 The stability of de Sitter for these models was also studied independently by A. Toporensky and col- 
laborators 1161 . 

■^As can be seen from the equations of motion, eqs. H16H and II17H are ill-defined in the \ — + limit. The 
case x = induces a relation between Ei and the other variables which makes it possible to eliminate S2 
and Si from the equations of motion. Therefore the case \ = has 4 fewer variables, which correspond to 
the last two pairs of eigenvalues given above. 
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We note that there are two zero eigenvalues. One of the zero eigenvalues just corresponds to 
the fact that this is a line of non-isolated equilibria. The second eigenvalue corresponds to a 
one-dimensional centre manifold and in order to determine the stability of the Kasner circle 
an analysis of this centre manifold is needed. Note also that the eigenvalue 2(1 — 2cos<^) 
arises from the curvature variable N. This is the same mode as in the GR case, where it 
causes vacuum type II transitions between points on the Kasner circle (see e.g. [T2"]). 
The centre manifold. Consider the invariant subspace given by (Q,T,, F,G) defined as 
follows: 

(£+,£_) = S(cos0,sin0), (S+i,S_i) = F£(cos </>, sin 0), 
(20) (£ + 2,£_ 2 ) = GS(cos^,sin0), N = B = fi A = 0, <\> = constant, 

and Q2 is given by the constraint equation. 
We consider the following perturbation: 

(Q, S, F, G) = (-3 + x, 1 + y, -3 + z, 18 + w). 

Close to the Kasner circle, the centre manifold can be parametrised by the linear combination 
defined by the variable 

and the equations of motion on the centre manifold turn into 

X' = X 2 + 0{X 3 ). 

This implies that close to X = 0, X is always increasing. In essence, the Kasner circle acts 
as a past attractor for some orbits, while for others it is unstable to the past. 

The centre manifold has a special importance for solutions asymptoting to the equilibrium 
point. Due to the zero-eigenvalue the decay of the solutions as they asymptote to the 
equilibrium point will be power-law (compared to exponential) in the dynamical time, r. 
Hence, at sufficiently early (or late) time, the decay rates will be dominated by the power-law 
behaviour. In the phase space this can be illustrated as follows. The solutions will approach 
the centre manifold exponentially rapidly, and then move along the centre manifold towards 
the equilibrium point in a power-law manner (see Fig[lJ. The centre manifold will therefore 
dominate the behaviour at sufficiently early (or late) times. 
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A comment on the results of ref. p~5] is in order. They also consider the stability of the 
Kasner solutions in quadratic theories of gravity, however, they conclude that the Kasner 
solutions are unstable. Technically, this is correct and they base their result on the presence 
of a logarithmic term in t. However, as we have seen here, there exist typical orbits having 
the Kasner solutions as a past attractor (in ref. [15j this amounts to choosing a negative k in 
their eq.(30) which would cause the singularity to appear for t — k > and thus hit would 
be finite there) . The cause of this problem is the centre manifold and it is our opinion that 
the dynamical system approach unravels the true nature of the (in)stability of the Kasner 
solutions more clearly. We therefore conclude that: in the space of solutions, there exists a 
set of non-zero measure which has the Kasner solutions as a past attractor. 

3.3. Quasi-FRW solution: X. This solution is given by 

E± = E±i = E ±2 =N = B = Q A = 0, Q = -2, Q 2 = 8. 

Its stability is determined by the eigenvalues 

l[x3], 2[x2], 3[x3], 4[x2]. 

The presence of all positive eigenvalues (ie unstable to the future) here implies that the 
solution is a past attractor. 

The metric corresponding to this equilibrium point is 

ds 2 = -dt 2 + t(dx 2 + dy 2 + dz 2 ), 

and has $ M „ — 0, as can be checked explicitly - it is the k — subcase of the radiation- like 
solution of eq. ^ and evolves like a flat radiation-dominated GR universe despite being a 
solution of the pure quadratic theory. 

3.4. Kasner-like solutions: JC(I). There is also a set of Kasner-like solutions when x < 
-1/2: 

Q = — 1, (E + , E_) = E(cos</>, sin<^), (E+i, E_i) = — £(cos</>, sin</>), 

(21) (E +2 ,E_ 2 ) = 2E(co S( />,sin0), £ 2 = l+^f* Q A = B = N = 0. 

-(1 + 2%) 

These equilibrium points have both both negative and positive eigenvalues and are therefore 
unstable. 

The metric corresponding to this exact solution of $ MJ/ = is 

(22) ds 2 = -dt 2 +t 2 \t- 4 ^dx 2 +t 2 ^++^-)dy 2 + t 2 ^-^ a -'>dz 2 

(a 2 + a 



-(1 + 2*)' 

3.5. Super-inflating FRW universe: E. This solution is given by 
Q = Q 2 = E± = E±i = E ±2 = N = fl A = B = 0. 
Its stability is determined by the eigenvalues: 




0[x2], -1, -3[x3], -- l±Wl+ p. 1 [x2] 



The exact solution corresponds to an inflating Einstein metric = Xg^v with A arbitrary; 
hence it is necessarily also a solution to = 0. 
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However, interestingly, this solution has a non-trivial centre manifold which causes a 
peculiar phenomenon. On the 2-dimcnsional centre manifold we can use B — X and Q\ = Y 
as coordinates. To lowest order in its neighbourhood we get 

X' « -X 2 , Y' « -XY. 
6 ' 6 

So if (4 — x)lx < 0, B < 0, this solution is stable, but is otherwise unstable. 

We can show that, by approximating the solution close to the equilibrium point, and 
assuming that the solution is stable, that the evolution of the scale factor is given by 
a(t) oc e H ° l I" 1 where t is the cosmological time. The Hubble scalar diverges linearly H = Hot 
and there are divergent curvature modes with R ~ t 2 as t — > oo. Hence, generic solutions 
with B < approaching E will have deceleration parameter q = — [1 + l/(Hot 2 )] < — 1 
with q — > — 1 at late times. These models are therefore 'super- inflationary' (or marginally 
'phantom-like'). 

3.6. Anisotropically-inflating type I universes: A(I). Unusually, for certain values of 
X and B, there are also exact solutions that describe anisotropic inflationary solutions of 
Bianchi type I: 

(£+,£_) = S(cos0,sin0), E 2 = „ 2 ( 4 ~ *) + B 
V + ' ' y ^' Vh 4(2 X +1) 

Q = E±i = S ±2 = N = 0. 
There are two classes of such solutions, depending on the values of B and £1a'- 

l : B = constant, il\ = — -, 

U A 8(2 X + 1)' 

(ii) : B = 0, Oa = constant. 

As long as x an d B take values for which S 2 > 0, these solutions exist. Their eigenvalues 
can be jointly expressed as: 



0[x2], -3[x3], -l(i±^Z*fj, -|(l±Vl 



1 + 8E 2 , -(1+4S+). 



Two of the zero eigenvalues appear because these are 2-parameter families of equilibrium 
points. However, since S 2 > 0, there will always be one unstable mode, and hence, these 
solutions are saddle points. 

The metrics corresponding to the case where B ^ can be written 

(23) ds 2 = -dt 2 + e 2bt [e-^+'dx 2 + e 2 (- + +^-)* d y 2 + e 2 ^"^-)*^ 2 

, 2 _ 1 + 8A(a +0) 2 2 _ 1 + 2A(4a + 0) 



6 2 = 



9/3 ' v + 18/3 



This set of solutions is defined so long as b 2 > and (<r? + a 2 _ ) > 0. As we can see, these 
solutions inflate anisotropically. Note that they are solutions of the theory in which A > 
but are not de Sitter spacetimes. Therefore, they show explicitly that the usual cosmic 
no-hair theorem of GR [T71 HS1 UHl HH1 HI] does not hold in these theories. In effect, the non 
linear ^ 0) terms contribute an effective stress tensor to the vacuum equations which 
violates the strong-energy condition needed for the cosmic no-hair theorem to hold in GR 
[22]. The essential features of this novel solution arise because of the contribution of the 
quadratic Ricci terms. If we put a = then there is no essential change in the character of 
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the solution but note that then solution does not have a GR limit when we take — > 0: it 
is an intrinsically quadratic solution of the higher-order theory. 

3.7. Anisotropically inflating type II universe: A(II). The phenomenon of anisotropic 
inflation is not confined to the Bianchi type I models. An anisotropically inflating solution 
of Bianchi type II is described by the critical point 

£+ = —-t, £- = £±i = £± 2 = Q = Q2 = 0, 
33 1 3 

(24) n A = ---N\ J B = 4(l + 8 X )iV 2 - i (ll-2 X ). 

This is the type II solution given in [9] . 

Its stability is determined by the eigenvalues 




-| (l±^a±Vb) , -3, -| (l± Vl + 16iV 2 

a = 15(3 - 167V 2 ), b = 9(9 - 624CW 2 - 243207V 4 ). 
The zero eigenvalue corresponds to the fact that this solution can be parametrised by B. 
Due to the presence of a positive eigenvalue, this solution is unstable to the future. 
There is a related set of equilibrium points for which 

£+ = ™, S_ = S ±1 = S ±2 = Q = Q2 = B = 0, 

(25) iV 2 = ^ ( VTqJ . fi A = constant. 

16 V 1 + 8% J 

These equilibrium points are defined for all Xj as long as N 2 > 0. 

3.8. Other solutions. We must stress that this list of equilibrium points of type II is not 
exhaustive. There are further critical points of the full dynamical system, which we do not 
consider here. 

4. Behaviours at infinity 

4.1. FRW case. Consider the behaviour at infinity for the FRW case; i.e., with £± = 
£±1 = E±2 = N — 0. We note that for large values of Q there are solutions which, to 
lowest order, give (noting that, after performing a translation of time, we can assume that 
Q diverges at r = without loss of generality) : 

(26) Q = ± Q 2 = ^#, b = ^, n A = ^f, 

It |t| 2 \t\ |t| 

where <52,0i Ah ^a.o ar e constants fulfilling the condition 

Boft\fi — —-. 

Hence, for these solutions at infinity to exist, the constants Bo and f^A.o need to have 
opposite sign; so, since > 0, we must have B < 0. 

Therefore, these behaviours at infinity come in two classes, according as to whether r < 0, 
which implies Q < and Q diverges into the future, or t > 0, which implies Q > and 
Q diverges to the past. These two choices for the sign of Q correspond to recollapsing and 
bouncing cosmological solutions, respectively. 
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Consider the case r < 0; by direct integration, we get for the Hubble scalar and the 
proper time t: 

H^H \r\K M'«^(t -i), 

and hence, H = (H 2 /2) (to — t) and the universe reaches a maximum size at to and contracts 
thereafter. The divergence of the variables for these cases is just a feature of the expansion- 
normalised variables. 

Let us now consider the FRW case in more detail. Assuming a FRW metric, the equations 
of motion can be written (without introducing expansion-normalised variables, but assuming 
(3a + /3) < 0): 

(27) 6H 2 H + 2HH - H 2 - r]H 2 + lu 2 = 0, 



where 



1 2 A 



? 7 = 777^ ; — 3T> u 



2(3a + /3)' 6(3a + /3)' 

and overdots denote d/dt. We note that u> can have either sign and rj > 0. 

We are interested in the case where there is a point where H = (i.e., the evolution has a 
turning point). Let us assume that this occurs for t = 0. Then we can find solutions which, 
close to t = 0, can be expanded as 

(28) H = ud+ ^yt 2 - |(6w - ?y)t 3 - |j(15w - r/)i 4 + 0(i 5 ). 

with H2 constant. As eq. (l2"T)) is ill-posed through H = 0, we have assumed that H(i) is 
analytic at t = 0. 

We note that in the case u> < 0, H goes from being positive to negative, hence the universe 
goes from an expanding to a collapsing phase. If to > 0, the universe goes from a collapsing 
to an expanding phase, and so the universe experiences a bounce at t = 0. We also note 
that this turning point is not symmetric with respect to t = as long as H2 ^ 0. The 
constant H2 is proportional to the constant Q2,o from above. Moreover, by comparing the 
approximate solution eq. (|28p . we can identify the behaviour at infinity, as given by eq. (|26p . 
as turning points of the expansion of the universe. Note that for the special case lu = 77 = 
the universe expands, reaches H = 0, and then continues to expand. 

4.2. General case. In the general Bianchi type I and II cases we can also find approximate 
solutions at infinity, with 

n - 1 n - ® 2 r - B ° o - ^ A0 
2t |r|5 \t\ \t\ 

( 29 ) L ± = 7-rr> L ±i = — r-j— , ^±2 = —-3-, A/ = — T . 

|rp |t| |r|a |t| 2 

The constants i?o, ^A0j S±,o, ^±1^0, S±2.o an d iVo must fulfill a complicated constraint 
arising from the 1/r 2 term of the constraint equation. 

From this behaviour, we can again show that these infinities correspond to turning points 
of the expansion of the universe between states of expansion and contraction. 

4.3. Transitions at infinity. Since the infinities described above correspond to turning 
points in the evolution, we can regularly pass through these infinities by switching to, for 
example, cosmological time. By choosing these infinities to occur at r = 0, we can pass 
through t — (t < h^t>0 will correspond to a transition from H > to H < 
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to the future, while the reverse will be a bounce in the past). In terms of the approxima- 
tions, eq. ()29|) we get the following transitions for the expansion- normalised variables and 
the dynamical time at infinity: 

(r, Q, Q 2 ,B, A , H±, E±i, E±2, N) ^ (-r, Q, -Q 2 , B, fi A> -E ± , E ±1) -E ±2 , -N), 

These transitions can now be used to continue these solutions through the formal infinity 
arising in state space. We should emphasise that for some invariant subspaces there might 
be different transitions to the one given above. 

5. Behaviour of FRW universes 

Let us consider the FRW case in further detail. It is defined by E± = E±i = E±2 = N = 
and we have f2 A > 0. It is also useful to define the following quantity: 

K=^, K' = 0. 

So if is a constant and we will subsequently replace all occurrences of J7 A with KB. Note 
also that B < implies K < 0, and B > implies K > 0. Moreover, we will use the 
constraint equation to solve for Qi. The equations of motion are then reduced to a two- 
dimensional system: 

(30) B' = -2QB, 

(31) Q' = -Iq2-3Q-Ib + ?-B 2 , 

which can now be analysed by using standard techniques. 

It is also useful to map the two-dimensional state space onto the compact two-dimensional 
unit disk, D 2 . This can be done by introducing polar coordinates (B, Q) — _R(cos0,sin</>), 
followed by, for example, the transformation: R i— > R/yl + R 2 . The behaviour of the state 
space at infinity is now mapped onto the unit circle S* 1 = 3D 2 . On the unit circle there are 
two points of special importance, labelled by the value of the angular variable, <j): 



Boo ■ tan < 




(32) Tloo : tan0 = +y^, <f> G (it, y 

The points B^ and IZoo correspond to the possible turning points of the expansion (bounce 
and recollapse, respectively). 

There are two other points on the unit circle that correspond to H — 0, namely 4> = ir/2 
and 4> = 37r/2, which will be denoted 7^+ and T^, respectively. Both of these points are 
in the invariant subspace B = 0, but are both unstable for models with B ^ 0. These two 
points correspond to u> = r\ = in eq. ([28|) . 

For B = 0, we can find the exact solution 

g(r) = Trk^> 

for which the metrics can be written 

ds 2 = fc2dr2 . + e 2T (dx 2 + dy 2 + dz 2 ) , e = -l,l. 

(l + ee -3r)| v ; = 

These metrics are solutions to $ M „ = and describe universes evolving as follows (see Fig(2]) : 
(1) e = 1: From J to E. 
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Figure 2. A sketch of the dynamical behaviour of vacuum FRW universes. 
For illustration, the dynamical system is mapped onto the unit disc as 
explained in the text. 

(2) e = — 1: For r < 0, from I to T M , while for r > 0, from T+ to E. Note that.after 
requiring analyticity, we have the transition T<^ — ► 7^+ at r = 0. 

For B = on the approach to E we have H = Ho at late times unlike the general case 
(B 0) for which the evolution gives H — Hot at late times. 

The behaviour of general FRW universes (B ^ 0) can be split into 5 different possibilities. 
The evolution can be as follows (see FigJ5]): 

(1) B > 

(2) B < 

(3) B < 

(4) B < 

(5) B < 

The last possibility suggests that it might even be possible for a universe to experience a 
sequence of expanding and collapsing phases (recall that a collapsing phase can be described 
by the same equations by running the dynamical time, r, backwards). 

6. Past behaviour of Bianchi-type universes 

In the analysis above we have seen that there are (at least) three possibilities for past 
behaviour of Bianchi models. Let us discuss each of these in turn. Again, we should 
emphasise that the list of equilibrium points for the Bianchi type II case is not exhaustive, 
and therefore behaviours other than these might occur in the general type II solution. 

6.1. Isotropic singularity. We have already seen that the radiation- like FRW universe 
denoted 2, is an attractor for the models studied here. This creates the possibility for an 
isotropic past singularity in these theories. Hence, the quadratic theories of gravity studied 



2 — 


dS. 


2 - 


E. 


2 - 




Boo 


-> E. 


Boo 
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here seem to have an isotropisation mechanism into the past provided by the quadratic 
curvature contributions. Let us give an argument that this isotropic singularity is indeed 
an attractor for all Bianchi models. 

Curvature modes. As explained earlier, the various Bianchi models can be described in terms 
of the structure coefficients via the symmetric tensor n a b and the vector a&. The evolution 
of these curvature variables are determined by the Jacobi identity; hence, the evolution of 
n a b and ab are therefore theory-independent. Introducing expansion-normalised variables in 
the standard manner, we obtain the equations of motion close to the isotropic singularity: 

Nl = -(Q + l)Ni, A' = -(Q + 1)A. 

Since Q = — 2 for the attractor I, we have 

N t m N i>0 e T , A « A e T 

and hence each quantity decays to the past r — > — oo. 

Shear modes. We do expect extra shear modes to appear; however, if we introduce the 
corresponding expansion-normalised variables for these modes, the evolution equations are, 
to linear order, the same as for E±. Consequently, this would give eigenvalues {1,2,3} for 
each of the additional shear modes. Therefore, we expect the shear to decay as 

Ei,- = Cije T + 0{e 2r ), 

at early times as r — > — oo. This is consistent with the situation found for the Bianchi type 
IX model by Cotsakis et al |15j . 

Perfect-fluid matter. With regards to the influence of matter on the evolution, let us assume 
the presence of a comoving perfect fluid with a barotropic equation of state: 

P = (7- l )P- 

Introducing an expansion-normalised energy-density, f2 = p/(3H 2 ), the conservation of 
energy yields the equation of motion: 

n' = -[2(Q + l) + (37-2)]ft. 

Hence, we have for the evolution of the fluid density with volume expansion, fl oc e' 4-37 ^ 1 ", 
and there is a change of behaviour of f2 at 7 = 4/3. However, in the evolution equations f2 
only appears as part of the product Bfl oc e^ 8 ~ 3j ^ T . This means that for all regular matter 
with 7 < 2, the matter term is subdominant and the isotropic singularity is stable to the 
past. In GR the situation is very different, and stability is only possible when 7 = 2 [24] . 

These conclusions have been drawn for the case of irrotational matter only. The situation 
with rotation or non-comoving 4- velocities, where u is not hypersurface orthogonal remains 
to be investigated. We know from experience with Bianchi-type universes in general relativ- 
ity that the introduction of all the possible non-comoving velocity components can produce 
a situation that is delicately sensitive to the equation of state because of the role played by 
the pressure and density in the conservation of angular momenta around the three orthogo- 
nal expansion axes, see refs [3SJH3]. The evolution of fluids possessing anisotropic pressures 
is also interesting and is likely to play a very important role in the late-time evolution of 
these universes, notably in the case where the anisotropic fluid has vanishing trace and is 
accompanied by an isotropic black-body radiation fluid. However, the way in which the 
quadratic Ricci terms have been found to mimic the effect of an isotropic radiation stress 
implies that the presence of a pure trace-free stress (for example a pure magnetic field) 
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may evolve in an influential way in these theories, just as was found in the case of general 
relativity [26| [27] . 

6.2. Kasner circle. The Kasner circle is an attractor for some orbits of type I and II and 
in these models this behaviour is therefore typical. The type II modes cause transitions 
between different points on the Kasner circle and since this is an exact solution for ordinary 
GR we expect chaos as we go to the more general vacuum type VIII, IX and Vl^/g models. 
However, even though chaotic solutions exist, we do not know whether this chaotic set is an 
attractor for typical orbits. 

6.3. Bounce. We have seen that there are bouncing solutions for the type I and II uni- 
verses. However, even though this behaviour seems to be typical for these models it is 
unclear whether they are typical for more general Bianchi models. The bouncing solutions 
correspond to solutions coming from infinity of the state space and are therefore particu- 
larly challenging to analyse. Apart from for the FRW case, rigorous results regarding the 
generality of these bouncing solutions are not known. The stability of the FRW behaviour 
in the type IX case [15] indicates that the closed FRW bounce behaviour will occur there 
but there may be other behaviours that are stable far from isotropy which have not yet 
been identified. The possibility for bounce is closely related to the possibility of recollapsing 
solutions in these theories. However, the analysis of this kind of behaviour using expansion- 
normalised variables is plagued by the fact that the state variables diverge. It does seem 
that recollapsing solutions are typical but a more rigorous analysis is lacking at the present 
time. 

7. Future behaviour 

The future behaviour of these cosmologies shows a dependence of the parameters of the 
theory. Unlike in general relativity, there does not exist a simple no-hair theorem when A > 
for the vacuum case, or where matter obeys the strong-energy condition. As discussed in 
our earlier paper, ref [9], this situation arises because the contribution of the stresses to 
the field equations ((2]) can mimic the form of a fluid with negative density and or pressures 
and so the energy-conditions assumed to hold in the general-relativistic no-hair theorems 
are effectively violated by the non-linear curvature terms. We showed that there exist 
vacuum solutions with A > which do not approach de Sitter at late times: they inflate 
anisotropically. A similar effect can be seen in the study of Kaloper into the effects of the 
Chern-Simons terms on the behaviour of a class of Bianchi type universes |28) . However, 
conditions can be identified for which de Sitter is an attractor at late times. 

7.1. De Sitter. For de Sitter to be an attractor we need B > and [l + 2A(4a + /3)]//3 < 0. 
However, even in this case, numerical simulations indicate that this behaviour is typical but 
not generic. But in the cases where B < 0, another late-time behaviour is possible: 

7.2. Super-inflating FRW. For B < and \ < or \ > 4 the equilibrium point E 
is a future attractor. Due to the presence of a non-trivial centre manifold, the universe 
generically approaches this point with a power-law evolution. As the solutions approach E, 
the universe eventually has H(t) = H t and q = -[1 + l/{H t 2 )}. 

Lately, there has been an interest for so-called "phantom cosmologies" with q < — 1 which 
can end in a singularity within finite time (a so-called 'big rip'). The evolution described 
above has eventually q < — 1 but approaches q = — 1 sufficiently fast to avoid the 'big 
rip' [29]. The deceleration parameter approaches —1 from below, while the Hubble scalar 
diverges linearly in cosmological time. Other types of finite-time singularity are possible in 
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which H and p remain finite but there is an infinity of the pressure and the acceleration 
[30] EH [32l [33] . These 'sudden' singularities occur without any violation of the strong- 
energy condition. It appears likely that they will occur in these quadratic theories also. The 
situation in the presence of lagrangians that involved only R was discussed in ref. [32] . 

8. Discussion 

In this paper we have provided a study of a range of anisotropic universes in a quadratic 
theory of gravity involving all allowed quadratic curvature invariants to appear in the la- 
grangian. We have seen that there is a possibility for a stable (or even generic) isotropic 
singularity in quadratic theories of gravity. There is a isotropic equilibrium point, repre- 
senting a FRW universe with scale factor a(t) oc y/t that is a solution of the vacuum field 
equations of the quadratic theory, which is stable into the past. Here, for simplicity, we have 
only considered Bianchi type I and type II models, but these incorporate both shear and 
3-curvature anisotropies, and an argument was given that this is also a past attractor for 
more general Bianchi models. Therefore, it appears that including quadratic terms provides 
us with a new mechanism for constraining the initial singularity to be isotropic. This is 
reminiscent of the situation proposed by Weyl curvature hypothesis, which envisages some 
measure of the Weyl curvature playing the role of a gravitational entropy, so it must initially 
be small (or zero, up to quantum corrections) in order to provide a subsequent gravitational 
arrow of time [34] [351 ES]. However, it must not be the case that the quadratic stresses 
drive the expansion towards isotropy on approach to any future singularity if the universe 
is closed and recollapses. 

This study raises many questions for further investigation of anisotropic and inhomoge- 
neous cosmologies that are more general than those considered here, which we will report 
elsewhere. We have also found that the presence of quadratic curvature terms leads to a 
range of new possibilities in situations which would be expected to yield simple de Sitter 
inflation in GR. The presence of a stress with p = —p need not produce future evolution to- 
wards the isotropic de Sitter metric locally. There is a possibility for inflationary behaviour 
which is anisotropic because of the effective stresses being contributed to the field equations 
by the quadratic terms. 

In ordinary GR, the general Bianchi types have been shown to possess a chaotic behaviour 
as we approach the initial singularity. This chaotic behaviour consists of a sequence of 
Kasner regimes connected via vacuum type II transitions and possibly frame rotations. On 
a dynamical time scale, the transitions are fairly short compared to the Kasner epochs 
[37l l38l [39] , Therefore, since we noticed that the Kasner circle does indeed act as a past 
attractor for some type I orbits, we could expect chaos to be present for typical orbits in the 
more general Bianchi models. Certainly, the chaotic vacuum type IX Mixmaster universe 
is a particular exact solution of the quadratic theory and its chaotic behaviour has been 
discussed in the presence of the quadratic terms [15] and it remains to be seen whether 
another form of chaotic evolution might arise, and be stable, far from the isotropic closed 
FRW radiation-like model. 
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